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We demonstrate that for gapped bilayer graphene, the nonlinear nature of the screening of an 
external disorder potential and the resulting inhomogeneity of the electron liquid are crucial for 
describing the electronic compressibility. In particular, traditional diagrammatic methods of many- 
body theory do not include this inhomogeneity and therefore fail to reproduce experimental data 
accurately, particularly at low carrier densities. In contrast, a direct calculation of the charge 
landscape via a numerical Thomas-Fermi energy functional method along with the appropriate bulk 
averaging procedure captures all the essential physics, including the interplay between the band gap 
and the inhomogeneity. 



I. INTRODUCTION 

Measurements of the electronic compressibility provide 
a way of characterizing the electron gas in both three- 
dimensional and two-dimensional materials and informa- 
tion about the nature of interactions between the elec- 
trons and the influence of the environment on the electron 
liquid can be gained. Therefore, it is highly important 
to have a clear theoretical understanding of experimental 
measurements of the compressibility. The compressibility 
K is given byi ^ ~ where n is the excess carrier 

density and fi is the chemical potential and so the key 
calculation is that of ^ . Recently, the compressibility of 
both monolayer^ and bilayer^ graphene has been exam- 
ined with capacitance probes, scanning single electron 
transistor (SET) microscopy and scanning tunneling 
microscopy (STM).^^^ g-,.^^ spatial resolu- 

tion of recent studies has shown a high degree of inho- 
mogeneity in the charge landscape of graphene systems, 
and revealed that the material used as the substrate 
has a significant impact on this inhomogeneity.^^ When 
the overall excess electronic density is close to zero (the 
so-called "charge- neutrality point"), the electron liquid 
breaks up into 'puddles' of electrons and holes, presum- 
ably to screen an external potential generated by disorder 
of some kind. Theoretical studies of graphene systems 
with charged impurities and corrugations or ripples 
have shown that either of these mechanisms may con- 
tribute to the observed inhomogeneity. 

Also, gapped electronic systems are highly important 
in many device applications, and bilayer graphene is an 
attractive material in this context since the band gap 
and carrier density can be controlled dynamically via 
gating^^i^ii^ Therefore a clear understanding of the inter- 
play between the inhomogeneity which is intrinsic to all 
graphene systems and the gapped nature of gated bilayer 
graphene is essential. In this paper we present a full anal- 
ysis of this issue via the theoretical consideration of the 
compressibility and comparison with recent experimen- 
tal work.^ Consistent with experimental findings (which 
are described in detail in Ref. [2I), we take the disorder 
to be arising from random quenched charged impurities 



in the environment of the bilayer graphene (BLG) with 
a two-dimensional (2D) impurity density of separated 
from the graphene layers by an average distance d. In 
Section HTH we apply the standard diagrammatic pertur- 
bation theory which is widely used to describe electron- 
impurity scattering in condensed matter systems. We 
show that this theoretical technique does not give the 
correct qualitative picture in the gapped regime when 
the inhomogeneity is strong. We contend that this fail- 
ure is due to the fact that this theory cannot incorporate 
the effects of the inhomogeneous charge distribution or 
the nonlinear nature of the screening. In Sec. IIVI we 
present a functional approach to calculate ^ based on 
and extending the Thomas-Fermi theory (TFT) of Ref. 
13, that is able to take into account the effect on the com- 
pressibility of the interplay of disorder and band-gap in 
the theoretically challenging regime when the band gap 
is of the same order, or smaller, than the strength of the 
disorder. 

We shall show that there are in fact two different cri- 
teria for assessing when the inhomogeneity is too strong 
for perturbative theories to be valid. The first is when 
the proportion of the graphene which is in the insulat- 
ing (and hence incompressible) state becomes significant. 
The second is when the average fluctuations character- 
ized by the root mean square of the density distribu- 
tion becomes large compared to the average carrier den- 
sity. This situation is qualitatively different in monolayer 
graphene^^^^^^ where there is no band gap and hence the 
screening nonlinearity has a much smaller effect on the 
compressibility since there is no mixed phase. 



II. THE CLEAN LIMIT 

In this section we give an overview of the single par- 
ticle physics of bilayer graphene in order to remind the 
reader of the most important points and to define our 
notation. The band structure of bilayer graphene can be 
approximated via a four-component Hamiltonian which 
describes the wave function amplitude on each of the four 
lattice sites in the unit cell.^ In this representation, there 
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are two branches in the conduction and valence bands, 
with one branch separated from the other by the inter- 
layer coupling energy 71 ^ 0.4eV. In the strongly inho- 
mogeneous regime, the split bands may become partially 
occupied even at low average carrier density and there- 
fore we keep all four bands in our analysis. The band 
structure is given by 



(1) 



where u is the band gap at/c = 0, i^ = +lin the conduc- 
tion band and = — 1 in the valence band, 6 = +1 in the 
split branches and h = —1 in the low-energy branches, vp 
is the Fermi velocity associated with monolayer graphene, 
and ?i = 1. These bands are each four- fold degener- 
ate due to the presence of spin and the two valleys in 
the Brillou in zone. The actual band gap^^ is given by 
u = uji I Y^^^t^^'^-O?. This band structure is illustrated in 
Fig. [TJa) for the ungapped case and two different band 
gaps. The quart ic (or 'sombrero') shape of the low-energy 
branches is clearly visible in the gapped examples. 



The compressibility associated with these bands can be calculated analytically by relating the density to the Fermi 
energy and computing the derivative.^ The full expression is 
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where VpTrn is the Fermi energy measured in terms of the density. 



Note that the first term implies that the electron liquid 
is incompressible at n = 0. There are three different 
cases because of the changes in the topology of the Fermi 
surface. At n the Fermi surface is ring shaped, but 
when the Fermi energy reaches the top of the sombrero 
part of the band structure (i.e. ja = u)^ this changes to a 
disk and hence there is a step in the value of Then, 
at much higher density (~ 2 x lO^^cm"^) the split band 
becomes occupied, and there are now two Fermi surfaces 
so that there is a second jump in The clean ^ is 

^ -"^ an dn 

shown in Fig. [l[b) for the three cases corresponding to 
the band structures in Fig. HJa). 



III. DIAGRAMMATIC PERTURBATION 
THEORY 
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FIG. 1. (a) The low-energy conduction and valence bands 
of bilayer graphene with and without a gate-induced band 
gap. (b) The clean given by Eq. ([2|). For both panels 
the dotted lines are = 0, the solid lines are u = lOOmeV, 
and the dashed lines are u = 200meV. The wave vector is 
measured in units of the lattice constant a. 



In this section we discuss perturbative methods for de- 
scribing the electron-impurity interaction. Inherent in 
this approach^^ is an average over disorder realizations 
which explicitely restores translational symmetry to the 
theory. We shall show that this approximation is not 
valid in the situation we discuss because it removes the 
possibility for inhomogeneity to form as a consequence 
of the electron-impurity scattering. In order to show 
that this failure is not an artifact of the specific level 
of approximation in the theory, we also apply the pertur- 
bation expansion keeping two different sets of diagrams, 
sum the infinite series associated with them, and find 



the same qualitative features in the predicted ^ which 
do not match the experimental results. 

In order to obtain ^ in this microscopic the- 
ory, the crucial feature which distinguishes the dis- 
ordered case from the clean case is the presence of 
the electron-impurity self-energy in the electron Green's 
function. We compute this self-energy within two differ- 
ent approximations^^— - the Born approximation (BA) 
and the self-consistent Born approximation (SCBA). In 
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the BA the self-energy is 

^ |y(k-k')pF„,,(k,k') 
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where Fc^Q,/(k, k') is the wave function overlap of the ini- 
tial and final states in the scattering process, Eo,\^ is the 
energy of an electron with wave vector k in band a from 
Eq. ([1]), 77 is a positive infinitesimal, and V^(k) is the 
screened impurity potential 
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with ^5 being the screening wave vector in the static ran- 
dom phase approximation, po the density of states of the 
clean system, and hi the dielectric constant. Note that 
the assumption of a homogeneous charge landscape also 
enters in the use of Qs as the screening wave vector. The 
SCBA takes into account the full Green's function for 
propagation between scattering events and for which the 
self-energy is given by 
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Note the self-consistent inclusion of the self-energy in the 
Green's function on the right-hand side. Once the self- 
energy has been obtained, the electron Green's function 
GaiKE) = [E- Ea\^ - J^aiKE)]'^ cau be straightfor- 
wardly computed. The density of states ^ is then re- 
lated to the imaginary part of the Green's function, since 



(6) 



where Qs = Qv = 2 are the spin and valley degeneracies, 
respectively. 

In Fig. [2] we show the calculated ^ for the clean case,^ 
the BA, and the SCBA as a function of the carrier den- 
sity. The right-hand panels show the same data on a 
logarithmic scale to emphasize the low-density features. 
In the ungapped case shown in Fig. Efa), the BA and 
SCBA give essentially the same result as the clean limit. 
When a band gap is present, as in Figs. [2fb) to[2fd), 
the clean limit shows a clear step occuring at the den- 
sity h = u'^/{v'^7r) which marks the density where the 
chemical potential leaves the sombrero region of the band 
structure and the topology of the Fermi surface changes 
from a ring to a disk.^ For n ^ n the BA and SCBA are 
similar to the clean limit, but for low-to- moderate den- 
sity n ^ h the strong modification of the density of states 
(DOS) near the band edge^^ implies that ^ is enhanced 
relative to the clean system but is still a decreasing func- 
tion as n becomes small. There is also a sharp divergence 
in the BA and SCBA for very small n which is not ob- 
served experimentallyj^i^ Additional structure for n < n 
in the BA comes from the non-trivial shape of the DOS 
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FIG. 2. (Color online) ^ for three different values of the band 
gap and impurity density for the diagrammatic perturbation 
theories. The points in panel (b) are taken from Fig. 2 of 
Ref. 4. Throughout, is given in units of eVnm^. 



near p = u in that approximation, which is smoothed 
out by the self consistency of the SCBA.^^ 

In Fig. Efb) we also show experimental data for the 
gapped regime. Note that the data in Ref. is a capaci- 
tance measurement, and in order to extract ^ we need 
accurate knowledge of experimental parameters, such as 
the impurity density, the gate-induced band gap, the 
stray capacitance, and the dielectric environment, all of 
which are known only approximately. Therefore the ex- 
perimental data shown here does not correspond to the 
parameters used in the calculation and hence we can- 
not expect quantitative agreement. Other experiments^!^ 
show the same qualitative features as in Ref. 4 although 
direct comparison to these data is not possible since the 
low-density ^ is obscured by the specifics of the experi- 
mental setup used in these measurements. We see in Fig. 
[2](b) that a broad peak forms in the experimental data 
at low density in complete qualitative contrast to the BA 
and SCBA theoretical results. Therefore, these theories 
utterly fail to capture the essential physics of the gapped 
system at low densities. This is, however, not unexpected 
since the low density regime is completely dominated by 
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the charged impurity induced random puddles of com- 
pressible and incompressible regions. 



IV. THOMAS-FERMI THEORY 

We now describe the TFT for the inhomogeneous sys- 
tem. In this approach, the carrier density landscape is 
obtained by minimizing a Thomas-Fermi energy func- 
tional of the spatially- varying density n(r) that includes a 
term due to the presence of disorder. The TFT is similar 
in spirit to the density functional theory (DFT)^ — but 
in TFT the kinetic energy operator is also replaced by 
a functional Exln]. This simplification makes the TFT 
valid only when the density profile varies on length scales 
larger than the Fermi wavelength, i.e., when |Vn/n| < 
/cf, where = \/7r|n| is the Fermi wave vector. This 
approach has been very successful in the context of trans- 
port calculations^ which provides strong phenomenolog- 
ical support for the use of this theory. Specifically for 
bilayer graphene, the puddle length scale is ~ 20nm and 
the density of carriers in the puddles is ~ lO^^cm"^ so 
that this inequality is marginally satisfied. However, we 
can also justify this approximation by pointing out that 
the root-mean-square of the density distribution is much 
larger than its average. At the charge- neutrality point 
(CNP) the average density n = cannot be taken as a 
measure of the typical carrier density inside the puddles 
and a better estimate is given by nrms- As a consequence, 
at low dopings (close to the CNP) nrms should be used in- 
stead of n in the inequality above. Given that nrms ~ 
we then conclude that the TFT is valid at all densities so 



long as Hi is not too small (n^ > lO^^cm"^). A full DFT 
for the disordered problem has been completed for mono- 
layer graphene^^ and shows very similar results to the 
TFT applied in the same context. However, the DFT 
is much more demanding of computational resources, and 
therefore it is not possible to simulate large lattice sizes 
or complete a comprehensive average over disorder real- 
izations in a reasonable timescale. Moreover, given the 
difficulty to quantitatively compare the theoretical and 
experimental results (since parameters, such as the im- 
purity density and stray capacitance, are not known ac- 
curately) and the complete failure of the diagrammatic 
methods to even achieve a gross qualitative description of 
the experimental results at low carrier density, our main 
motivation is to show that a functional method like the 
TFT is able to capture the qualitative features of the 
compressibility observed in experiments. The functional 
approach described below, notwithstanding the specifics 
of the TFT, is more than adequate to describe the com- 
pressibility of gapped systems in which the band gap is 
comparable to or smaller than the disorder strength. 
The TFT energy functional is given by 



E[n] = EK[n{r)] 
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where c^Vd / ^ is the bare disorder potential which is as- 
sumed to be due to the Coulomb interaction with ran- 
dom charged impurities with no spatial correlation and 
an equal probability of being positively or negatively 
charged. 



The first term is the kinetic energy where 
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is the ground state kinetic energy per excess carrier. 



The second term is the Hartree part of the electron- 
electron interaction, the third term is the contribution 
due to the disorder potential, and in the fourth term /i 
is the chemical potential. We neglect exchange and cor- 
relation terms^^ since, as we shall show below, ^ at low 
density is predominantly determined by the proportion 
of the sample which is incompressible and the inclusion of 
these terms will not change this. The ground state den- 
sity landscape is identified by the equation 6E/Sn = 0. 
Taking this variational derivative of E[n]^ we find 

^=e.[n(r)] + ^yrfr'^ + -Mr)-.. (9) 



Within our formalism it is fairly easy to assume the 
presence of spatial correlations among the impurities, and 
the presence of correlations has important effects on the 
transport properties of BLG;^^^^ however they do not 
modify the qualitative effects that the disorder has on 
the compressibility. 

Assuming that the clean ^ is valid locally, in Fig. [3] 
we show the spatial profile of the density distribution 
(left column) and ^ (right column) for a single realiza- 
tion of disorder with m = 3 x lO^^cm"^ for the gapless 
regime (first row) and the gapped case with u = 40meV 
and three values of the global charge density. The white 
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FIG. 3. (Color online) (Left column) Spatial density profile, 
and (Right column) spatial ^ profile for a single disorder 
realization. The units of the color bars are lO^^cm"^ for the 
density and eVnm^ for and = 3 x lO^^cm"^ through- 
out. The incompressible regions where n — are shown in 
white. The first row has u = and (n) = 0; the second row 
has u — 40meV and (n) — 0; the third row has u — 40meV 
and (n) ^ 2 X lO^^cm"^; and the fourth row has u = 40meV 
and {n) lO^^cm"^ 



regions are the parts of the graphene where the local 
density is zero and hence the graphene is incompress- 
ible. It is immediately noticeable that in the presence of 
a band gap (second row) there are large incompressible 
regions which are not present in the gapless case (even 
at zero excess carrier density), and these regions persist 
even when the average charge density (n) is significant 
((n) ^ 2 X lO^^cm"^, third row) and are still just visible 
when (n) ^ 4 x lO^^cm"^ (fourth row). More accurate 
functional methods than the TFT will not give substan- 
tially different values for the ratio of the sample that is 
covered by insulating (incompressible) regions, which we 
shall show is the dominant factor in determining ^ at 
low density. All these methods can do is to give slightly 
diflFerent values of n and density profiles inside metallic 




-4 -2 2 i 6 
n (10 cm ) 



(n) (lO^^cm"^) 



(c) 




<«>=lxl0'^ 



(d) 
0.151- 



0.05- 



1 


1 1 1 1 1 jy 1 








/ M=40,n=3xl0^^_ 




M=100, n~W^ 


^ 1 


M=200, R.= 10^^ 

1 1 1 1 1 1 1 



-4 -2 1^^ 6^ 8 10 
n (10 cm' ) 



"0 1 2 3 4 J 6.7 8 9 
{n) (10^ cm"^) 



FIG. 4. (Color online.) (a) P(n) at the charge neutrality 
point (n) = for different values of u and ra. For / 0, 
the arrow at n = represents the very narrow peak whose 
height is orders of magnitude larger than the y-scale used, 
(b) P(n — 0) for parameters corresponding to Fig. [2{b) to 
[2l^d). A dotted line at P(0) = 0.1 provided as a guide to 
the eye. (c) Evolution of P(n) with (n) for u — 40meV and 
Tii — ?>x lO^^cm"^. (d) /i as a function of doping. The legend 
in (d) also applies to panels (a) and (c). 



regions, both of which have very little infiuence on ^ for 
the situation of interest. 

By considering many disorder realizations we can cal- 
culate disorder averaged quantities and, in particular, the 
probability distribution function of the local carrier den- 
sity P{n). This function can then be used to compute 
the average density (n) = J n' P{n')dn' . As shown in Fig. 
m^a), P{n) is trimodal for (n) = in the gapped regime: 
it exhibits a large peak shown by an arrow at n = that 
quantifies the fraction of the sample ocupied by insulat- 
ing regions and two smaller and broader peaks centered 
around values of n which are determined by the nonlinear 
screening of the disorder potential and which therefore 
depend on both u and n^. Figure [D^ a) also shows that 
P(n) has a jump for n = h. As the doping increases, 
the fraction of the sample area covered by insulating re- 
gions P{n = 0) decreases as shown in Fig. IH^b). Notice 
the factor of 10^ difference in the vertical scale between 
Figs. Ilja) and [3Kb). The evolution of P(n) with (n) is 
shown in Fig. HJ^c). For finite doping, the distribution 
becomes asymmetric in n and becomes unimodal only at 
relatively large doping. In the unimodal regime, P(n) is 
approximated closely by a Gaussian centered around (n) . 

Experimental probes such as the capacitance mea- 
surements and scanning SET microscopy, simultaneously 
probe an area of the sample which is significantly larger 
than the puddles size as predicted in the TFT and mea- 
sured by STM.^^2 Therefore, an averaging procedure 
must be implemented to replicate ^ as a function of (n) 
as measured in those experiments. By disorder averaging 
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FIG. 5. (Color online.) (a) The predicted ^ in the TFT. (b) 
As in (a), but on a logarithmic scale to emphasize the low- 
density regime, (c) The density fluctuation due to disorder for 
the same parameters as in (a). The dotted line at rin^^s/^ — 
0.5 is included as a guide to the eye. 



the TFT results, we obtain the dependence of the average 
chemical potential (/i) (which is identical to fi in Eq. ([7])) 
with respect to the average density (n). Because of the 
non-linear screening, the relation between (/i) and (n) is 
also modified by the value of the gap and the strength 
of the disorder, as shown in Fig. ID^d). Thus, the TFT 
results clearly show the inhomogeneous nature of the car- 
rier density landscape in BLG in the presence of disorder. 

Using the TFT we calculate the average (^) = ^ 

which closely simulates the way in which (d/j^/dn) is ob- 
tained in both capacitance measurements^"— and in SET 
spectroscopy.^'^ 

Figure [5f a) shows the calculated ^ using the TFT for 
the same four sets of parameters as in Fig. [2l We imme- 
diately see that it exhibits qualitatively different behavior 
from the BA and SCBA, and at low density it shows a 
broad peak in qualitative agreement with the experimen- 
tal data. We stress that (as mentioned above) the data 
in Ref. is the only appropriate data for direct compar- 
ison with our theories, but since various parameters of 
the experimental system are unknown, we cannot expect 
quantitative agreement between our calculations and the 
measured data. Comparison of Fig. [2] with Figs. [5fa) 
and mb) shows that the deviation of the perturbation 
theory from the TFT occurs when P(0) > 0.1, indicating 



that the presence of insulating, incompressible regions is 
the dominating feature of the compressibility at low den- 
sity. In Fig. [5fc) we show our TFT calculated density 
fluctuation characterized by the root-mean-square value 
^rixis as a function of density. This clearly establishes 
that, when a band gap is present and as the fluctua- 
tions (i.e. inhomogeneity) become large with decreasing 
average density, the calculated TFT results for ^ start 
deviating substantially from the many-body perturbative 
ensemble averaged results, and when rirms/^ ~ 0.5, one 
must carry out the nonlinear screening theory to obtain 
the qualitatively correct features for the compressibility. 



V. DISCUSSION 

Therefore, we identify two different reasons for the fail- 
ure of the standard diagrammatic methods in gapped in- 
homogeneous systems. The first is the presence of strong 
inhomogeneity characterized by the parameter firms f^- 
When nrms/^ > cl with a ~ 1, the ground state cannot 
be assumed to be homogeneous and therefore implicit 
translational invariance incorporated in the disorder av- 
eraging step of the diagrammatic theory fails to quali- 
tatively describe the experimental situation. The exact 
value of a depends on the experimental quantity under 
consideration and the details of the experimental condi- 
tions. The second reason for failure is the existence of 
a random mixed inhomogeneous state where insulating 
(incompressible) and metallic (compressible) regions co- 
exist in the presence of a band gap of the order of or 
smaller than the disorder strength. The diagrammatic 
methods fail because they cannot account for this mixed 
state. For the compressibility, comparison of Fig. |4l^b) 
with Fig. [2] shows that the critical fraction of the sample 
area to be covered by insulating regions for the pertur- 
bative theories to give strong qualitative disagreement 
with the experiments is P(0) > 0.1. Thus, disorder has 
a much stronger qualitative effect at low carrier densities 
for gapped bilayer graphene^^^ than for the monolayer*^^ 
since BLG, by virtue of being a gapped system, can be 
in the random mixed state which is not accessible for 
ungapped systems. 

In conclusion, we have demonstrated that the nonlin- 
ear nature of the screening of an external disorder poten- 
tial in gapped bilayer graphene and the resulting charge 
inhomogeneity are crucial in understanding the ground 
state electronic properties for a wide range of experi- 
mentally relevant carrier densities. In particular, stan- 
dard many-body diagrammatic techniques assume that 
the density profile is homogeneous in both the screening 
and the Green's function, and therefore give qualitatively 
incorrect predictions for the compressibility in the pres- 
ence of an external band gap. In contrast, the TFT re- 
tains the inhomogeneity and non-linear screening of the 
density distribution in the energy functional E[n{T)] and 
therefore captures the essential physics of the system. 
We emphasize that this particular averaging procedure 
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discussed in Sec. |lVl simulating the experimental con- 
ditions, is simply inaccessible to any type of theories in- 
voking a homogeneous charge landscape to obtain many- 
body self-energy or broadening. We also point out that 
although we have presented results where the disorder 
potential is induced by random charged impurities, our 
general conclusion will remain valid for any form of disor- 
der which produces a scalar potential perturbation to the 
clean Hamiltonian, such as corrugations in the graphene 
sheet. Finally, although we have described calculation 
of the compressibility (and equivalently ^) the general 
logic of our argument applies to other observable quanti- 
ties also. For instance, application of the Kubo formula 
for transport with Green's functions derived in the same 
way as in Sec. IIIII will suffer from similar problems in 
the inhomogeneous regime. In this example, an effective 



medium theory^^ or a full quantum transport analysis 
that explicitly takes into account the inhomogeneities^ - 
should be used instead. 
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